Part 3: Making interactive maps with leaflet

Published

November 21, 2024

This is part 3 of the 3 part series.

Setup

Load the required R packages from CRAN and github.

From CRAN:

library(here)
library(data.table)
library(magrittr)
library(ggplot2)
library(leaflet)
library(sf)
library(devtools)
library(kableExtra)

From github:

devtools::install_github("deandevl/RspatialPkg")

Introduction

The benefit of creating a JavaScript map over a .jpg map as we did in our last post is that the map is “slippy,” that is, it slips around inside its container. You can drag to pan, scroll to zoom, click to show popups, etc. The downside, however, is that, since leaflet creates a JavaScript map, the map can only be shared in an interactive environment like a web browser. As such, leaflet is not a good choice for pasting images in papers and presentations, or for setting a snazzy new desktop background. Go back to Maps, Part 2 for that.

Review: Load data

The first dataset is a .geojson file containing geospatial descriptions of Philadelphia’s neighborhoods, courtesy of OpenDataPhilly. This dataset is polygon data and will form our basemap for layering on additional, more interesting, features.

The second dataset is geospatial point data, also provided by OpenDataPhilly, that contains information from the police department on shooting victims.

Read in the Philadelphia neighborhood data.

neighborhoods_file_path <- file.path(here(), "data", "philadelphia-neighborhoods.geojson")

neighborhoods_raw_sf <- sf::read_sf(neighborhoods_file_path)
sf::st_crs(neighborhoods_raw_sf) = 4326
Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -75.23049 ymin: 39.98491 xmax: -75.0156 ymax: 40.11269
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  NAME          LISTNAME MAPNAME Shape_Leng Shape_Area                  geometry
  <chr>         <chr>    <chr>        <dbl>      <dbl>        <MULTIPOLYGON [°]>
1 BRIDESBURG    Bridesb… Brides…     27815.  44586264. (((-75.06773 40.0054, -7…
2 BUSTLETON     Bustlet… Bustle…     48868. 114050424. (((-75.0156 40.09487, -7…
3 CEDARBROOK    Cedarbr… Cedarb…     20021.  24871745. (((-75.18848 40.07273, -…
4 CHESTNUT_HILL Chestnu… Chestn…     56394.  79664975. (((-75.21221 40.08604, -…
5 EAST_FALLS    East Fa… East F…     27401.  40576888. (((-75.18476 40.02829, -…
6 MOUNT_AIRY_E… Mount A… East M…     28846.  43152470. (((-75.18087 40.04325, -…

Read in the Philadelphia shootings data

shootings_file_path <- file.path(here(), "data", "philadelphia-shootings-2023.geojson")

shootings_raw_sf <- sf::read_sf(shootings_file_path)
sf::st_crs(shootings_raw_sf) = 4326
head(shootings_raw_sf)
Simple feature collection with 6 features and 22 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -75.15338 ymin: 40.02221 xmax: -75.12371 ymax: 40.04969
Geodetic CRS:  WGS 84
# A tibble: 6 × 23
  cartodb_id objectid  year dc_key   code  date_               time  race  sex  
       <int>    <int> <int> <chr>    <chr> <dttm>              <chr> <chr> <chr>
1          1 16074191  2023 2023350… 300   2023-05-28 20:00:00 04:1… B     M    
2          2 16074192  2023 2023350… 400   2023-06-02 20:00:00 17:0… B     M    
3          3 16074193  2023 2023350… 400   2023-08-05 20:00:00 23:1… B     M    
4          4 16074194  2023 2023350… 400   2023-06-03 20:00:00 22:3… B     M    
5          5 16074195  2023 2023350… 400   2023-08-05 20:00:00 23:1… B     M    
6          6 16074196  2023 2023350… 400   2023-08-05 20:00:00 23:1… W     M    
# ℹ 14 more variables: age <chr>, wound <chr>, officer_involved <chr>,
#   offender_injured <chr>, offender_deceased <chr>, location <chr>,
#   latino <int>, point_x <dbl>, point_y <dbl>, dist <chr>, inside <int>,
#   outside <int>, fatal <int>, geometry <POINT [°]>

Review: Clean data

leaflet requires that data be in WGS 84, so we would need to convert to WGS 84 (EPSG code: 4326) using sf::st_transform(shootings_raw_sf, crs = 4326) if it weren’t provided to us with that CRS.

If we want to run non-geospatial analysis on our shootings data, such as plotting shootings over time, calculating totals by demographic, and so on, we can drop the geospatial information and work with a standard tibble using sf::st_drop_geometry(shootings_raw_sf).

In the neighborhood data, rename a column.

neighborhoods_raw_sf <- data.table::as.data.table(neighborhoods_raw_sf) %>% 
  data.table::setnames(., old = "MAPNAME", new = "Label") %>% 
  sf::st_as_sf(.)

Geospatial layers in leaflet

You have the option of loading data either as the data = … argument in leaflet::leaflet() or waiting until subsequent layers to provide the data. As in our last post, we will add the data in each layer, since we are working with two distinct datasets.

Your first map

leaflet::leaflet() %>% 
  leaflet::addPolygons(data = neighborhoods_raw_sf) %>% 
  leaflet::addCircles(data = shootings_raw_sf)
Warning in validateCoords(lng, lat, funcName): Data contains 12 rows with
either missing or invalid lat/lon values and will be ignored

Basic map of Philadelphia gun violence using leaflet. Source: OpenDataPhilly.

Add a basemap

The leaflet package makes it easy to add map tiles, or “basemaps” to the layperson. You can either choose to call addTiles() with no arguments to get the default basemap from OpenStreetMap or choose to call addProviderTiles() to get one of the various third-party options. Our favorite is CartoDB.Voyager, but you can explore the entire set of options and pick your favorite.

leaflet::leaflet() %>% 
  leaflet::addProviderTiles(providers$CartoDB.Voyager) %>% 
  leaflet::addPolygons(data = neighborhoods_raw_sf) %>% 
  leaflet::addCircles(data = shootings_raw_sf)
Warning in validateCoords(lng, lat, funcName): Data contains 12 rows with
either missing or invalid lat/lon values and will be ignored

Leaflet map with provider tiles. Source: OpenDataPhilly.

Simple formatting adjustments

Let’s make some basic formatting adjustments to the polygons layer: line color, line weight, line opacity, and fill opacity (0 = no fill). We’ll also add a label, which will appear upon hover.

make_neighborhood <- function(){
  neighborhood_wig <- leaflet::leaflet() %>% 
    leaflet::addProviderTiles(providers$CartoDB.Voyager) %>% 
    leaflet::addPolygons(
      color = "#222",
      weight = 2,
      opacity = 1,
      fillOpacity = 0,
      label = ~lapply(Label, htmltools::htmlEscape),
      labelOptions = leaflet::labelOptions(direction = "top"),
      data = neighborhoods_raw_sf
    )
  return(neighborhood_wig)
}

make_neighborhood() %>% 
  leaflet::addCircles(data = shootings_raw_sf)
Warning in validateCoords(lng, lat, funcName): Data contains 12 rows with
either missing or invalid lat/lon values and will be ignored

Leaflet map adding a hover label, opacity, dark polygon borders.

Jitter points so that we can see them more clearly

You may have seen that some of the points in the plot above were darker than others. This is because we had overlapped multiple translucent circles. To avoid this issue, we will “jitter” our points, adding a small amount of random displacement in the x- and y-directions. To make this jitter consistent each time you render the plot, remember to set the seed value for the random jitter using set.seed().

set.seed(1776)
shootings_raw_sf <- sf::st_jitter(shootings_raw_sf, factor = 0.0004)

make_neighborhood() %>% 
  leaflet::addCircles(data = shootings_raw_sf)
Warning in validateCoords(lng, lat, funcName): Data contains 12 rows with
either missing or invalid lat/lon values and will be ignored

Leaflet map with jittered points.

Add labels for clearer communication

Our final set of aesthetic changes will be to our point layer. We add two new variables to our shootings dataset: a “Color” variable that encodes encodes the “fatal” variable into red and grey, and a “Popup” variable that summarizes key information about each shooting. This popup variable will appear in our map when we click on a point. In leaflet, labels appear upon hover, and popups appear upon click.

Create columns from shootings_raw_sf to define popup labels and point colors.

shootings_labels_sf <- data.table::as.data.table(shootings_raw_sf) %>% 
  .[, `:=`(
    Color = ifelse(fatal == 1, "#900", "#222"),
    Popup = paste0(
      "<b>", location, "</b>",
      "<br/><i>", date_, "</i>",
      "<br/><b>Race:</b> ", data.table::fcase(
        race == "B", "Black",
        race == "W", "White"
      ),
      "<br/><b>Sex:</b> ", data.table::fcase(
        sex == "M", "Male",
        sex == "F", "Female"
      ),
      "<br/><b>Age:</b> ", age,
      "<br/><b>Wound:</b> ", wound,
      "<br/><b>Fatal?:</b> ", data.table::fcase(
        fatal == 1, "Yes",
        fatal == 0, "No"
      )
    )
  )] %>% 
  .[, .(Color, Popup, geometry)] %>% 
  sf::st_as_sf(.)

Incorporate shootings_labels_sf over the neighborhood geometries.

make_neighborhood() %>% 
  leaflet::addCircles(
    color = ~Color,
    popup = ~Popup,
    data = shootings_labels_sf
  )
Warning in validateCoords(lng, lat, funcName): Data contains 12 rows with
either missing or invalid lat/lon values and will be ignored

Leaflet with shooting popup and point color

Choropleths in leaflet

Choropleths–maps in which each region is colored according to a summary statistic–are a powerful way to visualize data. In this example, let us suppose that we would like to show the total number of shootings in each neighborhood.

Join polygon data of neighborhoods_raw_sf with shootings_raw_sf.

neighborhoods_valid_raw_sf <- sf::st_make_valid(neighborhoods_raw_sf)
neigh_shootings_join_sf <- sf::st_join(neighborhoods_valid_raw_sf, shootings_raw_sf)

Convert the join to a data.table, group, and create “total_shootings” variable.

neighborhoods_raw_dt <- data.table::as.data.table(neighborhoods_raw_sf)

neigh_shootings_sf <- data.table::as.data.table(neigh_shootings_join_sf) %>% 
  .[, .(total_shootings = .N), by = Label] %>% 
  neighborhoods_raw_dt[., on = c("Label", "Label")] %>% 
  sf::st_as_sf(.)

Map the choropleth based on “total_shootings” variable.

pal <- leaflet::colorNumeric(
  "YlOrRd",
  domain = neigh_shootings_sf$total_shootings
)

leaflet::leaflet(neigh_shootings_sf) %>% 
  leaflet::addPolygons(
    color = "#222",
    weight = 2,
    opacity = 1,
    fillColor = ~pal(total_shootings),
    fillOpacity = 0.7
  )

Final map

In this final map, we add back our provider tiles, our label, and our highlight options, with no changes here from what had been done earlier in this post. We have also added a legend (and assigned the palette function to it), which describes the color range. Notice that we can change the opacity and location of the legend so that it is as unobtrusive as possible.

leaflet::leaflet(neigh_shootings_sf) %>% 
  leaflet::addPolygons(
    color = "#222",
    weight = 2,
    opacity = 1,
    fillColor = ~pal(total_shootings),
    fillOpacity = 0.7,
    label = ~lapply(Label, htmltools::HTML),
    labelOptions = leaflet::labelOptions(direction = "top"),
    highlightOptions = leaflet::highlightOptions(
      color = "#FFF",
      bringToFront = T
    )
  ) %>% 
  leaflet::addLegend(
    pal = pal,
    values = ~total_shootings,
    opacity = 0.7,
    title = "# shootings",
    position = "topleft"
  )

Leaflet final map with choropleth on shooting counts and neighborhood highlight

Conclusion

We would struggle to recreate and exact copy of ggplot2‘s maps in leaflet. But, that is to be expected. These two packages create two different types of maps–static and interactive–for different analytical purposes. What leaflet might lose in creating annotations and allowing for extremely precise aesthetic changes, it gains by allowing for paning, zooming, hovers, and popups.

Think of these two packages as complimentary tools in your analytics arsenal. Think carefully about when to use each one so that you can display data clearly, insightfully, and intuitively.